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A simple model of two-phase flow in porous media is presented. A connection is made to statistical 
mechanics by applying capillary power as a constraint. Stochastic sampling is then used to test the 
validity of this approach. Good agreement is found between stochastic sampling and time stepping 
for flow-rates above a transition value. 
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PACS numbers: 47.56.+r, 47.61. Jd, 89. 75. Da 

When one phase displaces another within a porous 
medium, complex patterns are known to emerge [l| . Less 
is known about the flow patterns formed under steady- 
state conditions. Dynamic effects, i.e., the dependence 
of the flow patterns on the total flow-rate also remain 
poorly understood. This is despite the significant im- 
portance to applications such as enhanced oil recovery, 
groundwater contamination and water transport in fuel 
cells [2-y|- Two- phase flow in porous media also holds 
theoretical interest as a complex system which exhibits 
self-organization. Self-organization is evidenced by capil- 
lary pressure drops adding up to reduce effective perme- 
ability rather than canceling out, even under steady-state 
conditions. 

Steady-state simulations have looked at the transport 
of disconnected oil ganglia 5] and relations between driv- 
ing pressure and fractional flow [6|], amongst other things. 
Experimental data can be classified by the model porous 
medium: networks etched in glass [7Hll|. Hele-Shaw cells 
12l Il3| and bead packings 14 1 . 

Two recent experiments have explored the relation ship 
between applied pressure drop AP and flow-rate |12l4l4 1 . 
They have reported a power law dependence, but have 
found different exponents. A two-dimensional network 
simulator has been used to explore the proposed power 
law |15l4l7| . In the simulator individual menisci are mod- 
eled; they are transported according to the flow field, and 
created and destroyed in a manner which crudely models 
snap-off and coalescence. Details are given in [lfj. 

This paper presents a model which captures similar 
behavior as the simulator without explicitly modeling 
menisci. The model consists of a network of capillaries 
and will be referred to as the Capillary Network Model 
(CNM). The capillaries model pore throats, which are 
the narrow connections between pores. A single variable 
f is assigned to each throat and a capillary pressure drop 
p c is given as a function of (p. 

The main control parameter used in simulations and 
experiments is the capillary number Ca = (iQ/a(j)A, 
where /1 is viscosity, Q is volumetric flow-rate, a is sur- 
face tension, (f> is porosity and A is the cross-sectional 
area of the sample. The effective permeability is the 
sum of the relative permeabilities, and may be defined as 
K c ff = Q/Qq where Qo is the flow-rate obtained by solv- 
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FIG. 1. In the CNM, each link in the network is a capil- 
lary containing a single droplet. The capillary models a pore 
throat; it is narrow in the middle and wide near the pores. 
The position of the droplet is given by tpi £ [0, 2-k]. 



ing Darcy's law for single-phase flow with pressure drop 
AP. K c g is found to be less than unity under steady- 
state conditions, which signifies that the mixture of two 
phases results in a larger resistance to flow than if only 
a sing le p hase is present. 

In [15|, [l6| a modified version of the Young-Laplace 
relation is used to obtain the capillary pressure drop of 
a single meniscus 
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where x is the position of the meniscus, £ is the length 
of the throat and r measures its width. In the CNM the 
explicit menisci are replaced by a single scalar variable for 
each pore throat. This variable ip is a coordinate which 
can take values between and 2n. A capillary pressure 
drop function is defined as 



p c (f) = — sin(7rs) sin tp, 



(2) 



where s is the non- wetting saturation of the throat, given 
by the length of the droplet divided by the length of the 
throat. For simplicity, all throats are assigned equal and 
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FIG. 2. (Color online) « e ff(Ca) for the simulator (black 
circles) and the CNM (red triangles). System size is 32x64. 



constant s. Eq. [2] is obtained by considering Eq. [T] for 
the case of two menisci forming a droplet within the pore 
throat, see Fig.Q] ip gives the position of the droplet. The 
size of the droplet determines the non- wetting saturation 
S, and can be removed as an independent parameter by 
redefining a. 

The porous medium is modeled as a 2D square lat- 
tice inclined at 45° relative to the main direction of flow. 
Boundary conditions are bi-periodic, such that the net- 
work may be mapped onto a torus. Flow is driven by a 
pressure drop applied across a cut through the network. 

Eq. [2] models pore throats that have a narrowing ge- 
ometry. In all other respects pore throats are considered 
to be cylindrical tubes with radius r. Network disorder 
may be introduced in r, I or both. In this work there is 
no disorder - all pore throats are equally wide and long. 
For all results given here, r = 0.2£ 

The Hagen-Poiseuille permeability for cylindrical tubes 
gives the flow-rate q by the Washburn equation, 

4 

q(Ap^) = -—(Ap + p c (ip)), (3) 

where \i is the effective viscosity of the phases contained 
in the throat and Ap is the pressure difference between 
the two pores connected by the throat. 

Eqs. [2] and [3] together with the network geometry and 
Kirchhoff's circuit laws define the CNM. 

Given a configuration consisting of {<fi} and Ca, the 
flow field is obtained by solving a system of linear equa- 
tions. The coordinates are then updated according to 
'■Pi = Qi/o. where a — irr 2 . Numerical time stepping is 
done by the Euler method, which is only first order accu- 
rate. The numerical error may be estimated by compar- 
ing the capillary power with the total dissipation. This 
ratio should be zero in the steady-state. It is found to be 
negligible for Ca > 10~ 3 , but grows to a few percent for 



Ca = 10 -4 . Below this the Euler method therefore does 
not provide reliable results. 

Fig. [2] shows simulation results for the CNM and a net- 
work simulator. In both cases Euler time stepping has 
been used. The network simulator is the same as that 
used in lallflj) but without disorder. Four independent 
runs underlie each data point. Initialization is done by 
a random initial configuration and a gradual increase in 
surface tension (gradually decreasing Ca). Saturation is 
a trivial parameter in the CNM, but not in the simula- 
tor. The simulator results are for S = 0.4. To facilitate 
comparison between the CNM and the simulator, s in 
Eq. [5] was set equal to 0.4. The two models can be seen 
to produce similar, but not identical results. 

Consider Eq. [3] The right-hand side consists of two 
terms, one proportional to Ap and one to p c . In obvious 
notation q = qo + q c - Multiplying with q and divid- 
ing by —irr 4 /8£fj, gives a statement of conservation and 
conversion of energy, d = do + d c , where d is the heat 
dissipated through viscous shear, always negative by def- 
inition, —do = —Ap q is equivalent to the power provided 
by a pump driving a flow q with an external pressure drop 
Ap and d c — p c q is capillary dissipation. Replacing —do 
with w, where w represents the energy of the pump, and 
—d c with capillary power w c gives w + w c + d = 0. This 
states that the absolute value of heat dissipation within 
a single throat is equal to the sum of applied power and 
capillary power. 

For the porous medium as a whole the capillary power 
adds up to W c = ^ u> c . The pump power is W = —APQ 
and the total heat dissipation is D = J2 d- In general, 
the pump power and the heat dissipation are not equal. 
The difference must be due to the capillaries, so we have 
W c = —W — D for the total capillary power. If W c is 
positive the capillary power adds to the pump power to 
increase heat dissipation; it becomes an energy source. If 
W c is negative the capillary power absorbs some of the 
pump power to decrease heat dissipation; it becomes an 
energy sink. 

Fig. El shows the development of W, —D and W c dur- 
ing a single simulation run. Initially, for a completely 
random configuration, W c is positive: the capillaries re- 
lease energy. W c approaches zero as steady-state is ap- 
proached, after which both pump power and dissipation 
fluctuate around the same average value. 

A thermodynamics of two-phase flow in porous media 
was first suggested in [18]. There, total dissipation was 
suggested as being analogous to energy. The preceeding 
discussion implies that instead of total dissipation it is 
capillary power which is constrained and thus provides a 
connection to statistical mechanics. 

In the following, tools from statistical mechanics usu- 
ally reserved for classical equilibrium will be applied. The 
underlying idea is to consider the steady-state as gov- 
erned by a balance between drive and dissipation - a 
dissipative equilibrium. (W c ) = states that, on aver- 
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FIG. 3. (Color online) Time development of W (black, 
initially increasing), — D (red, initially decreasing) and W c 
(inset) for an initially random configuration. The x-axis of the 
inset is a zoom-in compared to the main figure. Ca = 10 -3 is 
held constant and the system size is 128x256. 
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FIG. 4. (Color online) A(Ca) for the CNM (plusses). In- 
set: 1 — K c ff (black triangles) and k c h (red squares) with the 
same x-axis as the main figure. The data are from the same 
simulations as the CNM results from Fig. [2] 



age, pump power must equal heat dissipation, which is a 
requirement for the steady-state. 

Constructing a space of eigenstates {fi}, where I in- 
dexes eigenstates, gives the number of microscopic con- 
figurations for a state {a;} as 



C({ai}) 



N\ 



ai!a 2 !a 3 ! . 



(4) 



where a/ is the occupancy number of eigenstate I and 
^2 ai = N is the number of coordinates [19( . In light of 
the developments so far, and in order to establish the gov- 
erning principle of steady-state two-phase flow in porous 
media, the capillary power W c is applied as a constraint 
on Eq.|31 Using a Lagrange multiplier A for the constraint 
gives the partition function 



* = £• 



-xje 



(5) 



where the sum runs over the eigenstates of (p. J%? is 

^ = -£><* (6) 

i 

where the sum runs over the number of coordinates, p c = 
p c {(fl) and q is 

q(ifl) = dq q p{q\tpi). (7) 



The pdf 's p may be obtained from time stepping at a 
given Ca. Capillary power 



W C (X) = -Z- x d x Z = -NZ^J^PcV e Xp ^ 



(8) 



is a strictly monotonous function of A, with only a single 
solution A for a given value of W c . Fig. |4] gives A(Ca) 
obtained by solving Eq.[8]for W c = 0. A transition occurs 
at Ca « 0.1. This value is in the followingreferred to as 
Ca up , to emphasize the connection with 16 1. 

The inset of Fig. |4] shows that above Ca up , (1 — n c s) 
scales with Ca. The best fit to a power law gives (1 — 
«cff) ~ Ca~ 2 ' 06 . A further discussion of scaling laws 
above Ca up is given in connection with the mean-field 
argument. 

It should be noted that W C (X) = implies that d\Z = 
0, i.e., dissipative equilibrium implies that the partition 
function has an extremal value with respect to variations 
of the constraint. For W C (X) — to be valid, the con- 
straint must apply equally for each configuration within 
the ensemble. Temporal correlations in the fluctuations 
of W c are ignored by this approach. 

Using a Metropolis algorithm 20j, [2l| and the value of 
A obtained from time stepping it is possible to generate 
an ensemble of configurations governed by the constraint 
of W c = 0. The resulting ensemble will be dominated 
by the most probable state (according to Eq. [4j which 
satisfies the constraint. A comparison with the ensemble 
obtained by time stepping constitutes a non-trivial test 
of the validity of the statistical mechanics approach. 

Stochastic sampling is done by replacing some ran- 
domly chosen coordinates with new, random coordinates. 
After a trial update, AW C is calculated as the change in 
W c after the update. The new configuration is accepted 
if AIU C is negative, and with probability e ~ AAW/c oth- 
erwise. This is a standard Metropolis algorithm, where 
W c is analogous to energy and A is analogous to inverse 



temperature. 

Fig. [5]shows K ff(Ca) and 1 — K e ff(Ca) from both time 
stepping and stochastic sampling. For Ca above Ca up 
the results are in agreement. Below Ca up , time stepping 
and stochastic sampling produce different results. 

Inspection of time stepping simulations suggest that 
above Ca up , q may be approximated by 



<?(</?) =Qo -Qc sin (p, 



(9) 



where go is a mean flow-rate and g c is the maximum per- 
turbation caused by a capillary pressure drop. At Ca up , 
9c ~ 9o- Below Ca up time stepping simulations show 
that the dependence of q on ip is no longer sinusoidal. 
Eq. [§] represents a mean-field solution: a homogeneous 
flow field with perturbations that only depend on the 
local variable. 

Applying Eq. [9] to Eq. [5] allows an analytic deter- 
mination of the scaling of Ca and (I — k c s) with A. 
To obtain this, the Boltzmann factor is expanded as 
e ~Ajsr K ^ _ \j>ff : omitting terms of A 2 and higher or- 
der. Terms with odd powers of sinip sum to zero. From 
Ca ~ (g), Ca ~ go is obtained. Considering (w c ) = 0, 
Ca ~ A -1 / 2 results. Finally, k cS ~ {q} 2 /{q 2 } gives 
(1 — n e e) ~ A. This gives (I — « c ff) ~ Ca -2 . 

From time stepping and above Ca up , the CNM gives 
Ca ~ A" 048 , (1 - n cS ) ~ A " and (1 - k cS ) ~ Ca" 206 , 
see the inset of Fig. [5] and Fig. El 

In summary, a simple model of two-phase flow in 
porous media has been presented. It has been found to 
produce results for effective permeability that are quali- 
tatively similar to a more detailed simulator. Capillary 
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FIG. 5. (Color online) Comparison of K e ff (Ca) obtained from 
time stepping (black squares) and stochastic sampling (blue 
plusses). Inset: 1 — K e « with the same x-axis as the main 
figure. The solid red line is a power law with exponent —2.06. 
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FIG. 6. (Color online) Ca(A). The solid red line is a power 
law with exponent —0.48. Inset: 1 — K e g with the same x-axis 
as the main figure. The solid red line is a power law with 
exponent 0.99. 



power has been defined and used as a constraint to pro- 
duce a partition function. At high flow-rates, time step- 
ping and stochastic sampling have been shown to pro- 
duce the same ensemble, given a constraint of zero cap- 
illary power. Scaling exponents obtained from a mean- 
field theory are in agreement with the numerical results 
at these high flow-rates. At lower flow-rates, stochastic 
sampling does not reproduce the results from time step- 
ping. 

The main result of this work is to identify zero capil- 
lary power as the constraint which governs steady-state 
two-phase flow in porous media. This is nothing but a 
statement of conservation of energy. 

Discussions with S. Sinha and G. T0ra are gratefully 
acknowledged. 
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